Prior Analysis
Susannah Tysor
20/03/2020
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(ggplot2)
library(bayesplot)## This is bayesplot version 1.7.1## - Online documentation and vignettes at mc-stan.org/bayesplot## - bayesplot theme set to bayesplot::theme_default()## * Does _not_ affect other ggplot2 plots## * See ?bayesplot_theme_set for details on theme settinglibrary(rstan)## Loading required package: StanHeaders## rstan (Version 2.19.2, GitRev: 2e1f913d3ca3)## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)library(cowplot)##
## ********************************************************## Note: As of version 1.0.0, cowplot does not change the## default ggplot2 theme anymore. To recover the previous## behavior, execute:
## theme_set(theme_cowplot())## ********************************************************library(testthat)##
## Attaching package: 'testthat'## The following object is masked from 'package:dplyr':
##
## matcheslibrary(parallel)
library(purrr)##
## Attaching package: 'purrr'## The following object is masked from 'package:testthat':
##
## is_nullrstan_options(auto_write=TRUE)# generate samples from a truncated normal distribution. n = how many samples, mean and sd are mean and sd, and min and max are the limits of the distribution/truncation points. possibly my first ever use of while loops. Don't @ me.
rtnorm <- function(n, mean, sd, min, max) {
x <- rnorm(n, mean=mean, sd=sd)
x <- x[x >= min & x <= max]
while(length(x) < n) {
newx <- rnorm(1, mean=mean, sd=sd)
while(newx <= min | newx >= max) {
newx <- rnorm(1, mean=mean, sd=sd)
}
x <- c(x, newx)
}
length(x)==n
return(x)
}
# simulate data from an ordinal logistic model and format it as input for stan. input is a list for stan, as is output.
simulate_data <- function(input) {
simu <- rstan::stan(file='dirichlet prior/covar_sim.stan', iter=1, chains=1, algorithm="Fixed_param", data=input)
simu_params <- rstan::extract(simu)
input_data_for_model <- list("N" = input$N, "K" = input$K, "x" = input$x, "y" = array(simu_params$y[1,]))
return(input_data_for_model)
}
# plot simulated data for a sanity check.
simplot <- function(datalist) {
inputdf <- data.frame(datalist)
p1 <- ggplot(inputdf, aes(x=x, y=y)) +
geom_jitter(shape=1, height=0.2) +
geom_vline(xintercept = h) +
ggtitle("Simulated data with cutpoints")
print(p1)
p2 <- ggplot(inputdf, aes(x=x, colour=as.factor(y))) +
stat_ecdf() +
geom_vline(xintercept=h) +
theme(legend.position = "none") +
ggtitle("Cumulative x for 3 states")
print(p2)
p3 <- cowplot::plot_grid(p1,p2)
print(p3)
}
# append the correct shape (gamma prior on cutpoints) and rate (exponential prior on beta) parameters to the simulated data list and fit a model with a gamma prior. simdat is a list of datasets simulated from stan models and pars is a vector of additional parameters. Relies on global params beta_slow and beta_fast.
fit_gamma_model <- function(simdatlist, pars) {
#choose whether to use data simulated with a rapid or slow transition
if (pars$beta == beta_slow) {
simdat <- simdatlist$simdat_slow
}
if (pars$beta == beta_fast) {
simdat <- simdatlist$simdat_fast
}
#extract parameters for prior distribtuions
simdat$shape <- pars$shape
simdat$rate <- pars$beta_rate
#fit the model
fitgam <- stan(file='dirichlet prior/gamma/gamma_covar.stan', data=simdat, chains=4)
return(fitgam)
}
## append a label (string) to all columnnames in a dataframe (x)
label_names <- function(x, label) {
colnames(x) <- paste0(colnames(x), "_", label)
return(x)
}
# function to plot modeled parameters with label being a string to label the graph (usually prior and whether groups were included and pardf is the modeled parameter dataframe) trueparams is a one row dataframe of true parameters. h is a global parameter have fun!
parplot <- function(pars) {
cutplot <- mcmc_areas(pars, regex_pars="c.\\d_model") +
geom_vline(xintercept=c(pars$c.1_true, pars$c.2_true))
betaplot <- mcmc_areas(pars, pars="beta_model") +
geom_vline(xintercept=pars$beta_true)
h1plot <- mcmc_areas(pars, regex_pars = "h1") +
geom_vline(xintercept=h[1])
h2plot <- mcmc_areas(pars, regex_pars = "h2") +
geom_vline(xintercept=h[2])
allplot <- cowplot::plot_grid(cutplot, betaplot, h1plot, h2plot,
nrow=1, ncol=4,
labels=paste("model", unique(pars$modelid_true),
"rate=", unique(pars$beta_rate_true),
"shape=", unique(pars$shape_true)))
print(allplot)
}
# function to calculate the difference between modeled parameters (in the model_params dataframe) and true parameters (h is globally declared). Where true params is a dataframe
posterior_differencer <- function(pars) {
c1_diff <- pars$c.1_model - pars$c.1_true
c2_diff <- pars$c.2_model - pars$c.1_true
h1_diff <- pars$h1_model - h[1]
h2_diff <- pars$h2_model - h[2]
beta_diff <- pars$beta_model - pars$beta_true
diffframe <- data.frame(c1_diff, c2_diff, h1_diff, h2_diff, beta_diff)
return(diffframe)
}
# function to plot histograms of differences between true params and modeled params (model_params dataframe).
diffplotter <- function(diffs, pars) {
#diffs <- posterior_differencer(pars)
cuts <- mcmc_areas(diffs, regex_pars = "c") +
ggtitle("", subtitle = "differences between modeled and true params")+
xlim(c(-15,15))
opars <- mcmc_areas(diffs, pars=c("h1_diff", "h2_diff", "beta_diff")) +
xlim(c(-1,1))
cowplot::plot_grid(cuts, opars, labels=paste("model", unique(pars$modelid_true),
"rate=", unique(pars$beta_rate_true),
"shape=", unique(pars$shape_true)))
}Prior choice
In my original phenology model, I use
- a gamma distribution for the prior on the cutpoints
- an exponential distribution distribution for the transition speed/slope \(\beta\)
- normal distributions for group effects on the slope
Michael Betancourt thinks that an induced dirichlet prior may be a good choice for ordinal logistic models.
Normally, priors are bottom up - you choose a distribution for each prior. This is tricky for the cutpoints in an ordinal logistic model because they’re defined on an abstract latent space so you can’t easily use your domain expertise, but a good prior is incredibly important in these models, because ordered logistic models with covariates are inherently non-identifiable. (Because \(beta\)s and cutpoints depend on one another.)
Betancourt says that > To avoid the non-identifiability of the interior cut points from propagating to the posterior distribution we need a principled prior model that can exploit domain expertise to consistently regularize all of the internal cut points at the same time. Regularization of the cut points, however, is subtle given their ordering constraint. Moreover domain expertise is awkward to apply on the abstract latent space where the internal cut points are defined.
Goals
- Understand the dynamics of the gamma and induced dirichlet priors and 2) compare them
Simulate data
First let’s simulate some data that’s kind of like mine.
We have three potential states K and a latent effect/covariate x. Data is most likely to be collected around state 2. x is always positive.
I’ll simulate 2 datasets - one with a fast transition speed (\(\beta = 2\)) and one with a slow transition speed (\(\beta = 0.5\)). I want the transitions to occur at the same x for both datasets. The halfway transition points h will be at x= 8 and x=12. So the cutpoints c for the slow transition will be at 4 and 6 and for the fast transition at 16 and 24.
Simulate datasets
N <- 500
K <- 3
beta_slow <- 0.5
beta_fast <- 2
cutpoints_slow <- c(4,6)
cutpoints_fast <- c(16,24)
h <- cutpoints_slow/beta_slow #half transition points, engineered to be identical for slow and fast transitions
x <- rtnorm(n=N, mean=mean(h), sd=2, min=0, max=20) #covariate
testthat::test_that("half transition points identical", {
testthat::expect_equal(cutpoints_slow/beta_slow, cutpoints_fast/beta_fast)
})inputs_for_sim_slow <- list("N" = N, "K" = K, "c" = cutpoints_slow, "beta"=beta_slow, "h" = h, "x" = x)
inputs_for_sim_fast <- list("N" = N, "K" = K, "c" = cutpoints_fast, "beta"=beta_fast, "h" = h, "x" = x)
# simulate data, graph it, and prepare data for model fitting
simdat_slow <- simulate_data(input=inputs_for_sim_slow)##
## SAMPLING FOR MODEL 'covar_sim' NOW (CHAIN 1).
## Chain 1: Iteration: 1 / 1 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0 seconds (Warm-up)
## Chain 1: 0.000153 seconds (Sampling)
## Chain 1: 0.000153 seconds (Total)
## Chain 1:slowsimplot <- simplot(simdat_slow)
simdat_fast <- simulate_data(input=inputs_for_sim_fast)##
## SAMPLING FOR MODEL 'covar_sim' NOW (CHAIN 1).
## Chain 1: Iteration: 1 / 1 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0 seconds (Warm-up)
## Chain 1: 0.00015 seconds (Sampling)
## Chain 1: 0.00015 seconds (Total)
## Chain 1:fastsimplot <- simplot(simdat_fast)
cowplot::plot_grid(slowsimplot, fastsimplot, nrow=2, labels=c("slow", "fast"), vjust=2)
Recapture parameters: Gamma
Now I want to try to recapture parameters.
I need to know how good I have to be at choosing the shape parameter for the gamma prior on the cutpoints to recapture parameters as well. So I’ll choose a prior for gamma that roughly centers it between the cutpoints, which I think is a “good” option (though given the long right tail, maybe a lower value would be better). Then halve and double it.
I also want to understand how the exponential prior on beta affects the ability to recapture parameters. I’ll try rates of 1,2,3.
Gamma prior with covariate
beta_rate <- c(1:3) # rate parameters for exponential prior on beta
# shape parameters for gamma prior on cutpoints, for slow and fast transitions
shape_slow <- mean(cutpoints_slow)
shape_slow_small <- shape_slow/2
shape_slow_big <- shape_slow*2
slowshapes <- c(shape_slow, shape_slow_small, shape_slow_big)
shape_fast <- mean(cutpoints_fast)
shape_fast_small <- shape_fast/2
shape_fast_big <- shape_fast*2
fastshapes <- c(shape_fast, shape_fast_big, shape_fast_small)
slowframe <- data.frame(beta = beta_slow, shape = slowshapes, beta_rate=beta_rate, c.1 = cutpoints_slow[1], c.2=cutpoints_slow[2], h1=h[1], h2=h[2])
fastframe <- data.frame(beta=beta_fast, shape= fastshapes, beta_rate=beta_rate, c.1 = cutpoints_fast[1], c.2=cutpoints_fast[2], h1=h[1], h2=h[2])
parframe <- rbind(slowframe, fastframe) %>%
tidyr::expand(tidyr::nesting(beta, c.1, c.2, shape), h1, h2, beta_rate)
parframe$modelid <- 1:nrow(parframe) #label the models
parframe <- dplyr::select(parframe, modelid, beta, c.1, c.2, h1, h2, shape, beta_rate)
knitr::kable(parframe, caption="beta parameters used to simulate data for a slow (0.5) and fast (2) transion. And shape and rate parameters to be used in the priors on beta and cutpoints to attempt to recover parameters")| modelid | beta | c.1 | c.2 | h1 | h2 | shape | beta_rate |
|---|---|---|---|---|---|---|---|
| 1 | 0.5 | 4 | 6 | 8 | 12 | 2.5 | 1 |
| 2 | 0.5 | 4 | 6 | 8 | 12 | 2.5 | 2 |
| 3 | 0.5 | 4 | 6 | 8 | 12 | 2.5 | 3 |
| 4 | 0.5 | 4 | 6 | 8 | 12 | 5.0 | 1 |
| 5 | 0.5 | 4 | 6 | 8 | 12 | 5.0 | 2 |
| 6 | 0.5 | 4 | 6 | 8 | 12 | 5.0 | 3 |
| 7 | 0.5 | 4 | 6 | 8 | 12 | 10.0 | 1 |
| 8 | 0.5 | 4 | 6 | 8 | 12 | 10.0 | 2 |
| 9 | 0.5 | 4 | 6 | 8 | 12 | 10.0 | 3 |
| 10 | 2.0 | 16 | 24 | 8 | 12 | 10.0 | 1 |
| 11 | 2.0 | 16 | 24 | 8 | 12 | 10.0 | 2 |
| 12 | 2.0 | 16 | 24 | 8 | 12 | 10.0 | 3 |
| 13 | 2.0 | 16 | 24 | 8 | 12 | 20.0 | 1 |
| 14 | 2.0 | 16 | 24 | 8 | 12 | 20.0 | 2 |
| 15 | 2.0 | 16 | 24 | 8 | 12 | 20.0 | 3 |
| 16 | 2.0 | 16 | 24 | 8 | 12 | 40.0 | 1 |
| 17 | 2.0 | 16 | 24 | 8 | 12 | 40.0 | 2 |
| 18 | 2.0 | 16 | 24 | 8 | 12 | 40.0 | 3 |
So there are 2 simulated datasets (beta=0.5 or beta=2) and I’m going to try to fit them both with 3 rates for the beta’s exponential prior and 3 shapes for the cutpoints’ gamma prior, a total of 18 model runs.
parlist <- split(parframe, seq(nrow(parframe))) # format parframe so it works with parLapply better
names(parlist) <- paste0("m", 1:nrow(parframe))
datlist <- list(simdat_slow=simdat_slow, simdat_fast=simdat_fast) # list of simulated datasets
# run all models, parallelized
# make a cluster, leaving 20 cores free for other folks
no_cores <- parallel::detectCores() - 20
cl <- parallel::makeCluster(no_cores)
# export the stuff you need to run on the cluster
parallel::clusterExport(cl, c("fit_gamma_model", "parlist", "beta_fast", "beta_slow", "datlist"))
parallel::clusterEvalQ(cl, c(library(rstan), library(StanHeaders)))## [[1]]
## [1] "rstan" "ggplot2" "StanHeaders" "stats" "graphics"
## [6] "grDevices" "utils" "datasets" "methods" "base"
## [11] "rstan" "ggplot2" "StanHeaders" "stats" "graphics"
## [16] "grDevices" "utils" "datasets" "methods" "base"
##
## [[2]]
## [1] "rstan" "ggplot2" "StanHeaders" "stats" "graphics"
## [6] "grDevices" "utils" "datasets" "methods" "base"
## [11] "rstan" "ggplot2" "StanHeaders" "stats" "graphics"
## [16] "grDevices" "utils" "datasets" "methods" "base"
##
## [[3]]
## [1] "rstan" "ggplot2" "StanHeaders" "stats" "graphics"
## [6] "grDevices" "utils" "datasets" "methods" "base"
## [11] "rstan" "ggplot2" "StanHeaders" "stats" "graphics"
## [16] "grDevices" "utils" "datasets" "methods" "base"
##
## [[4]]
## [1] "rstan" "ggplot2" "StanHeaders" "stats" "graphics"
## [6] "grDevices" "utils" "datasets" "methods" "base"
## [11] "rstan" "ggplot2" "StanHeaders" "stats" "graphics"
## [16] "grDevices" "utils" "datasets" "methods" "base"
##
## [[5]]
## [1] "rstan" "ggplot2" "StanHeaders" "stats" "graphics"
## [6] "grDevices" "utils" "datasets" "methods" "base"
## [11] "rstan" "ggplot2" "StanHeaders" "stats" "graphics"
## [16] "grDevices" "utils" "datasets" "methods" "base"
##
## [[6]]
## [1] "rstan" "ggplot2" "StanHeaders" "stats" "graphics"
## [6] "grDevices" "utils" "datasets" "methods" "base"
## [11] "rstan" "ggplot2" "StanHeaders" "stats" "graphics"
## [16] "grDevices" "utils" "datasets" "methods" "base"
##
## [[7]]
## [1] "rstan" "ggplot2" "StanHeaders" "stats" "graphics"
## [6] "grDevices" "utils" "datasets" "methods" "base"
## [11] "rstan" "ggplot2" "StanHeaders" "stats" "graphics"
## [16] "grDevices" "utils" "datasets" "methods" "base"
##
## [[8]]
## [1] "rstan" "ggplot2" "StanHeaders" "stats" "graphics"
## [6] "grDevices" "utils" "datasets" "methods" "base"
## [11] "rstan" "ggplot2" "StanHeaders" "stats" "graphics"
## [16] "grDevices" "utils" "datasets" "methods" "base"
##
## [[9]]
## [1] "rstan" "ggplot2" "StanHeaders" "stats" "graphics"
## [6] "grDevices" "utils" "datasets" "methods" "base"
## [11] "rstan" "ggplot2" "StanHeaders" "stats" "graphics"
## [16] "grDevices" "utils" "datasets" "methods" "base"
##
## [[10]]
## [1] "rstan" "ggplot2" "StanHeaders" "stats" "graphics"
## [6] "grDevices" "utils" "datasets" "methods" "base"
## [11] "rstan" "ggplot2" "StanHeaders" "stats" "graphics"
## [16] "grDevices" "utils" "datasets" "methods" "base"
##
## [[11]]
## [1] "rstan" "ggplot2" "StanHeaders" "stats" "graphics"
## [6] "grDevices" "utils" "datasets" "methods" "base"
## [11] "rstan" "ggplot2" "StanHeaders" "stats" "graphics"
## [16] "grDevices" "utils" "datasets" "methods" "base"
##
## [[12]]
## [1] "rstan" "ggplot2" "StanHeaders" "stats" "graphics"
## [6] "grDevices" "utils" "datasets" "methods" "base"
## [11] "rstan" "ggplot2" "StanHeaders" "stats" "graphics"
## [16] "grDevices" "utils" "datasets" "methods" "base"
##
## [[13]]
## [1] "rstan" "ggplot2" "StanHeaders" "stats" "graphics"
## [6] "grDevices" "utils" "datasets" "methods" "base"
## [11] "rstan" "ggplot2" "StanHeaders" "stats" "graphics"
## [16] "grDevices" "utils" "datasets" "methods" "base"
##
## [[14]]
## [1] "rstan" "ggplot2" "StanHeaders" "stats" "graphics"
## [6] "grDevices" "utils" "datasets" "methods" "base"
## [11] "rstan" "ggplot2" "StanHeaders" "stats" "graphics"
## [16] "grDevices" "utils" "datasets" "methods" "base"
##
## [[15]]
## [1] "rstan" "ggplot2" "StanHeaders" "stats" "graphics"
## [6] "grDevices" "utils" "datasets" "methods" "base"
## [11] "rstan" "ggplot2" "StanHeaders" "stats" "graphics"
## [16] "grDevices" "utils" "datasets" "methods" "base"
##
## [[16]]
## [1] "rstan" "ggplot2" "StanHeaders" "stats" "graphics"
## [6] "grDevices" "utils" "datasets" "methods" "base"
## [11] "rstan" "ggplot2" "StanHeaders" "stats" "graphics"
## [16] "grDevices" "utils" "datasets" "methods" "base"
##
## [[17]]
## [1] "rstan" "ggplot2" "StanHeaders" "stats" "graphics"
## [6] "grDevices" "utils" "datasets" "methods" "base"
## [11] "rstan" "ggplot2" "StanHeaders" "stats" "graphics"
## [16] "grDevices" "utils" "datasets" "methods" "base"
##
## [[18]]
## [1] "rstan" "ggplot2" "StanHeaders" "stats" "graphics"
## [6] "grDevices" "utils" "datasets" "methods" "base"
## [11] "rstan" "ggplot2" "StanHeaders" "stats" "graphics"
## [16] "grDevices" "utils" "datasets" "methods" "base"
##
## [[19]]
## [1] "rstan" "ggplot2" "StanHeaders" "stats" "graphics"
## [6] "grDevices" "utils" "datasets" "methods" "base"
## [11] "rstan" "ggplot2" "StanHeaders" "stats" "graphics"
## [16] "grDevices" "utils" "datasets" "methods" "base"
##
## [[20]]
## [1] "rstan" "ggplot2" "StanHeaders" "stats" "graphics"
## [6] "grDevices" "utils" "datasets" "methods" "base"
## [11] "rstan" "ggplot2" "StanHeaders" "stats" "graphics"
## [16] "grDevices" "utils" "datasets" "methods" "base"# fit the models
gammafits <- parallel::parLapply(cl, parlist, function(x) {fit_gamma_model(simdatlist = datlist, pars=x)})
parallel::stopCluster(cl) #close the cluster
# extract params
paramsgam <- lapply(gammafits,
function(x) {data.frame(rstan::extract(x) ) } )
#param summary
sum <- lapply(paramsgam, summary)
#bind true and model pars even tho it will make a giant df
paramsgam <- map(paramsgam, label_names, label="model")
parlist <- map(parlist, label_names, label="true")
paramsgam <- map2(paramsgam, parlist, cbind)
#plot params and diffs
map(paramsgam, parplot)
## $m1
##
## $m2
##
## $m3
##
## $m4
##
## $m5
##
## $m6
##
## $m7
##
## $m8
##
## $m9
##
## $m10
##
## $m11
##
## $m12
##
## $m13
##
## $m14
##
## $m15
##
## $m16
##
## $m17
##
## $m18
map(paramsgam, posterior_differencer) %>%
map2(.y=paramsgam, .f=diffplotter)## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.## Warning: Removed 2 rows containing missing values (geom_segment).## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.## Warning: Removed 2 rows containing missing values (geom_segment).## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.## Warning: Removed 2 rows containing missing values (geom_segment).## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.## Warning: Removed 2 rows containing missing values (geom_segment).## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.## Warning: Removed 2 rows containing missing values (geom_segment).## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.## Warning: Removed 2 rows containing missing values (geom_segment).## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.## Warning: Removed 2 rows containing missing values (geom_segment).## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.## Warning: Removed 2 rows containing missing values (geom_segment).## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.## Warning: Removed 1 rows containing missing values (geom_segment).## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.## $m1
##
## $m2
##
## $m3
##
## $m4
##
## $m5
##
## $m6
##
## $m7
##
## $m8
##
## $m9
##
## $m10
##
## $m11
##
## $m12
##
## $m13
##
## $m14
##
## $m15
##
## $m16
##
## $m17
##
## $m18
# calculate whether true value is in 50% hpdi
# kludge HPDI
HPDIlow <- function(x, prob) {
HPDI <- rethinking::HPDI(x, prob=prob)
return(HPDI[1])
}
HPDIhigh <- function(x, prob) {
HPDI <- rethinking::HPDI(x, prob=prob)
return(HPDI[2])
}
calc_HPDI <- function(params, prob) {
low <- params %>% dplyr::summarise_at(vars(ends_with("model")), HPDIlow, prob=0.5)
high <- params %>% dplyr::summarise_at(vars(ends_with("model")), HPDIhigh, prob=0.5)
# awkward formatting
hdpis <- dplyr::full_join(low, high) %>%
select(-contains("lp"))
colnames(hdpis) <- stringr::str_replace(colnames(hdpis), "_model", "")
# true param
true <- params %>% dplyr::summarise_at(vars(ends_with("true")), unique)
colnames(true) <- stringr::str_replace(colnames(true), "_true", "")
true <- select(true, colnames(hdpis))
# more awkward formatting
compframe <- dplyr::full_join(hdpis, true) %>%
t(.) %>%
data.frame()
colnames(compframe) <- c("low", "high", "true")
compframe$params <- rownames(compframe)
# true param in interval?
tf <- compframe %>% mutate(inint = true > low & true < high)
return(tf)
}
in50 <- map(paramsgam, calc_HPDI, prob=0.5)## Joining, by = c("c.1_model", "c.2_model", "beta_model", "h1_model", "h2_model",
## "lp___model")## Warning: Column `c.1_model` has different attributes on LHS and RHS of join## Warning: Column `c.2_model` has different attributes on LHS and RHS of join## Warning: Column `beta_model` has different attributes on LHS and RHS of join## Warning: Column `h1_model` has different attributes on LHS and RHS of join## Warning: Column `h2_model` has different attributes on LHS and RHS of join## Warning: Column `lp___model` has different attributes on LHS and RHS of join## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")## Joining, by = c("c.1_model", "c.2_model", "beta_model", "h1_model", "h2_model",
## "lp___model")## Warning: Column `c.1_model` has different attributes on LHS and RHS of join## Warning: Column `c.2_model` has different attributes on LHS and RHS of join## Warning: Column `beta_model` has different attributes on LHS and RHS of join## Warning: Column `h1_model` has different attributes on LHS and RHS of join## Warning: Column `h2_model` has different attributes on LHS and RHS of join## Warning: Column `lp___model` has different attributes on LHS and RHS of join## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")## Warning: Column `c.1_model` has different attributes on LHS and RHS of join## Warning: Column `c.2_model` has different attributes on LHS and RHS of join## Warning: Column `beta_model` has different attributes on LHS and RHS of join## Warning: Column `h1_model` has different attributes on LHS and RHS of join## Warning: Column `h2_model` has different attributes on LHS and RHS of join## Warning: Column `lp___model` has different attributes on LHS and RHS of join## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")## Warning: Column `c.1_model` has different attributes on LHS and RHS of join## Warning: Column `c.2_model` has different attributes on LHS and RHS of join## Warning: Column `beta_model` has different attributes on LHS and RHS of join## Warning: Column `h1_model` has different attributes on LHS and RHS of join## Warning: Column `h2_model` has different attributes on LHS and RHS of join## Warning: Column `lp___model` has different attributes on LHS and RHS of join## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")## Warning: Column `c.1_model` has different attributes on LHS and RHS of join## Warning: Column `c.2_model` has different attributes on LHS and RHS of join## Warning: Column `beta_model` has different attributes on LHS and RHS of join## Warning: Column `h1_model` has different attributes on LHS and RHS of join## Warning: Column `h2_model` has different attributes on LHS and RHS of join## Warning: Column `lp___model` has different attributes on LHS and RHS of join## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")## Warning: Column `c.1_model` has different attributes on LHS and RHS of join## Warning: Column `c.2_model` has different attributes on LHS and RHS of join## Warning: Column `beta_model` has different attributes on LHS and RHS of join## Warning: Column `h1_model` has different attributes on LHS and RHS of join## Warning: Column `h2_model` has different attributes on LHS and RHS of join## Warning: Column `lp___model` has different attributes on LHS and RHS of join## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")## Warning: Column `c.1_model` has different attributes on LHS and RHS of join## Warning: Column `c.2_model` has different attributes on LHS and RHS of join## Warning: Column `beta_model` has different attributes on LHS and RHS of join## Warning: Column `h1_model` has different attributes on LHS and RHS of join## Warning: Column `h2_model` has different attributes on LHS and RHS of join## Warning: Column `lp___model` has different attributes on LHS and RHS of join## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")## Warning: Column `c.1_model` has different attributes on LHS and RHS of join## Warning: Column `c.2_model` has different attributes on LHS and RHS of join## Warning: Column `beta_model` has different attributes on LHS and RHS of join## Warning: Column `h1_model` has different attributes on LHS and RHS of join## Warning: Column `h2_model` has different attributes on LHS and RHS of join## Warning: Column `lp___model` has different attributes on LHS and RHS of join## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")## Warning: Column `c.1_model` has different attributes on LHS and RHS of join## Warning: Column `c.2_model` has different attributes on LHS and RHS of join## Warning: Column `beta_model` has different attributes on LHS and RHS of join## Warning: Column `h1_model` has different attributes on LHS and RHS of join## Warning: Column `h2_model` has different attributes on LHS and RHS of join## Warning: Column `lp___model` has different attributes on LHS and RHS of join## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")## Warning: Column `c.1_model` has different attributes on LHS and RHS of join## Warning: Column `c.2_model` has different attributes on LHS and RHS of join## Warning: Column `beta_model` has different attributes on LHS and RHS of join## Warning: Column `h1_model` has different attributes on LHS and RHS of join## Warning: Column `h2_model` has different attributes on LHS and RHS of join## Warning: Column `lp___model` has different attributes on LHS and RHS of join## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")## Warning: Column `c.1_model` has different attributes on LHS and RHS of join## Warning: Column `c.2_model` has different attributes on LHS and RHS of join## Warning: Column `beta_model` has different attributes on LHS and RHS of join## Warning: Column `h1_model` has different attributes on LHS and RHS of join## Warning: Column `h2_model` has different attributes on LHS and RHS of join## Warning: Column `lp___model` has different attributes on LHS and RHS of join## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")## Warning: Column `c.1_model` has different attributes on LHS and RHS of join## Warning: Column `c.2_model` has different attributes on LHS and RHS of join## Warning: Column `beta_model` has different attributes on LHS and RHS of join## Warning: Column `h1_model` has different attributes on LHS and RHS of join## Warning: Column `h2_model` has different attributes on LHS and RHS of join## Warning: Column `lp___model` has different attributes on LHS and RHS of join## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")## Warning: Column `c.1_model` has different attributes on LHS and RHS of join## Warning: Column `c.2_model` has different attributes on LHS and RHS of join## Warning: Column `beta_model` has different attributes on LHS and RHS of join## Warning: Column `h1_model` has different attributes on LHS and RHS of join## Warning: Column `h2_model` has different attributes on LHS and RHS of join## Warning: Column `lp___model` has different attributes on LHS and RHS of join## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")## Warning: Column `c.1_model` has different attributes on LHS and RHS of join## Warning: Column `c.2_model` has different attributes on LHS and RHS of join## Warning: Column `beta_model` has different attributes on LHS and RHS of join## Warning: Column `h1_model` has different attributes on LHS and RHS of join## Warning: Column `h2_model` has different attributes on LHS and RHS of join## Warning: Column `lp___model` has different attributes on LHS and RHS of join## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")## Warning: Column `c.1_model` has different attributes on LHS and RHS of join## Warning: Column `c.2_model` has different attributes on LHS and RHS of join## Warning: Column `beta_model` has different attributes on LHS and RHS of join## Warning: Column `h1_model` has different attributes on LHS and RHS of join## Warning: Column `h2_model` has different attributes on LHS and RHS of join## Warning: Column `lp___model` has different attributes on LHS and RHS of join## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")## Warning: Column `c.1_model` has different attributes on LHS and RHS of join## Warning: Column `c.2_model` has different attributes on LHS and RHS of join## Warning: Column `beta_model` has different attributes on LHS and RHS of join## Warning: Column `h1_model` has different attributes on LHS and RHS of join## Warning: Column `h2_model` has different attributes on LHS and RHS of join## Warning: Column `lp___model` has different attributes on LHS and RHS of join## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")## Warning: Column `c.1_model` has different attributes on LHS and RHS of join## Warning: Column `c.2_model` has different attributes on LHS and RHS of join## Warning: Column `beta_model` has different attributes on LHS and RHS of join## Warning: Column `h1_model` has different attributes on LHS and RHS of join## Warning: Column `h2_model` has different attributes on LHS and RHS of join## Warning: Column `lp___model` has different attributes on LHS and RHS of join## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")## Warning: Column `c.1_model` has different attributes on LHS and RHS of join## Warning: Column `c.2_model` has different attributes on LHS and RHS of join## Warning: Column `beta_model` has different attributes on LHS and RHS of join## Warning: Column `h1_model` has different attributes on LHS and RHS of join## Warning: Column `h2_model` has different attributes on LHS and RHS of join## Warning: Column `lp___model` has different attributes on LHS and RHS of join## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")print(in50)## $m1
## low high true params inint
## 1 3.2501236 3.9272662 4.0 c.1 FALSE
## 2 4.9079239 5.6276046 6.0 c.2 FALSE
## 3 0.4066015 0.4748422 0.5 beta FALSE
## 4 7.8773879 8.2708234 8.0 h1 TRUE
## 5 11.7804592 12.1964284 12.0 h2 TRUE
##
## $m2
## low high true params inint
## 1 3.1747532 3.8290328 4.0 c.1 FALSE
## 2 4.9277306 5.6227433 6.0 c.2 FALSE
## 3 0.4031203 0.4692065 0.5 beta FALSE
## 4 7.8930745 8.2767062 8.0 h1 TRUE
## 5 11.8244379 12.2359328 12.0 h2 TRUE
##
## $m3
## low high true params inint
## 1 3.1904716 3.8115469 4.0 c.1 FALSE
## 2 4.9024314 5.5857277 6.0 c.2 FALSE
## 3 0.4086681 0.4723668 0.5 beta FALSE
## 4 7.8722504 8.2621884 8.0 h1 TRUE
## 5 11.8230640 12.2186223 12.0 h2 TRUE
##
## $m4
## low high true params inint
## 1 3.4753182 4.1142259 4.0 c.1 TRUE
## 2 5.2798238 5.9367518 6.0 c.2 FALSE
## 3 0.4363857 0.4991233 0.5 beta FALSE
## 4 8.0089769 8.3662952 8.0 h1 FALSE
## 5 11.7793695 12.1459087 12.0 h2 TRUE
##
## $m5
## low high true params inint
## 1 3.4433367 4.0945892 4.0 c.1 TRUE
## 2 5.1325837 5.8360014 6.0 c.2 FALSE
## 3 0.4246286 0.4904984 0.5 beta FALSE
## 4 8.0189219 8.3836955 8.0 h1 FALSE
## 5 11.7668734 12.1422131 12.0 h2 TRUE
##
## $m6
## low high true params inint
## 1 3.370690 4.0211574 4.0 c.1 TRUE
## 2 5.186302 5.8966790 6.0 c.2 FALSE
## 3 0.418465 0.4850303 0.5 beta FALSE
## 4 7.999701 8.3563466 8.0 h1 TRUE
## 5 11.791000 12.1658031 12.0 h2 TRUE
##
## $m7
## low high true params inint
## 1 3.9334486 4.5523476 4.0 c.1 TRUE
## 2 5.7414360 6.4048148 6.0 c.2 TRUE
## 3 0.4744432 0.5376042 0.5 beta TRUE
## 4 8.2128315 8.5227315 8.0 h1 FALSE
## 5 11.6305653 11.9513433 12.0 h2 FALSE
##
## $m8
## low high true params inint
## 1 3.9378608 4.557596 4.0 c.1 TRUE
## 2 5.7729105 6.447308 6.0 c.2 TRUE
## 3 0.4702384 0.533441 0.5 beta TRUE
## 4 8.2315496 8.543833 8.0 h1 FALSE
## 5 11.6784857 11.993660 12.0 h2 FALSE
##
## $m9
## low high true params inint
## 1 3.9731254 4.5753073 4.0 c.1 TRUE
## 2 5.6989545 6.3508008 6.0 c.2 TRUE
## 3 0.4805663 0.5419001 0.5 beta TRUE
## 4 8.1834663 8.5009195 8.0 h1 FALSE
## 5 11.6774571 11.9984272 12.0 h2 FALSE
##
## $m10
## low high true params inint
## 1 11.159384 12.494009 16 c.1 FALSE
## 2 18.014987 19.826380 24 c.2 FALSE
## 3 1.448803 1.602843 2 beta FALSE
## 4 7.672810 7.841917 8 h1 FALSE
## 5 12.192592 12.356752 12 h2 FALSE
##
## $m11
## low high true params inint
## 1 11.185234 12.402412 16 c.1 FALSE
## 2 17.852498 19.520479 24 c.2 FALSE
## 3 1.445014 1.587621 2 beta FALSE
## 4 7.656411 7.817452 8 h1 FALSE
## 5 12.170768 12.333943 12 h2 FALSE
##
## $m12
## low high true params inint
## 1 10.848001 12.049589 16 c.1 FALSE
## 2 17.692949 19.370861 24 c.2 FALSE
## 3 1.421928 1.564056 2 beta FALSE
## 4 7.664332 7.829884 8 h1 FALSE
## 5 12.197229 12.360644 12 h2 FALSE
##
## $m13
## low high true params inint
## 1 12.524896 13.899434 16 c.1 FALSE
## 2 19.839749 21.731258 24 c.2 FALSE
## 3 1.604186 1.765372 2 beta FALSE
## 4 7.760157 7.913378 8 h1 FALSE
## 5 12.146768 12.297673 12 h2 FALSE
##
## $m14
## low high true params inint
## 1 12.320277 13.634744 16 c.1 FALSE
## 2 19.447600 21.265906 24 c.2 FALSE
## 3 1.592442 1.745851 2 beta FALSE
## 4 7.735040 7.883800 8 h1 FALSE
## 5 12.148636 12.299303 12 h2 FALSE
##
## $m15
## low high true params inint
## 1 12.213129 13.531782 16 c.1 FALSE
## 2 19.267255 21.137336 24 c.2 FALSE
## 3 1.575907 1.733495 2 beta FALSE
## 4 7.729892 7.877785 8 h1 FALSE
## 5 12.143976 12.292382 12 h2 FALSE
##
## $m16
## low high true params inint
## 1 15.299240 16.740585 16 c.1 TRUE
## 2 23.408511 25.483513 24 c.2 TRUE
## 3 1.930559 2.106695 2 beta TRUE
## 4 7.854474 7.987955 8 h1 FALSE
## 5 12.089570 12.224809 12 h2 FALSE
##
## $m17
## low high true params inint
## 1 15.322033 16.795263 16 c.1 TRUE
## 2 23.497265 25.605068 24 c.2 TRUE
## 3 1.939864 2.115915 2 beta TRUE
## 4 7.845753 7.980946 8 h1 FALSE
## 5 12.070039 12.206998 12 h2 FALSE
##
## $m18
## low high true params inint
## 1 14.975978 16.394557 16 c.1 TRUE
## 2 23.018115 24.980390 24 c.2 TRUE
## 3 1.907683 2.074526 2 beta TRUE
## 4 7.856500 7.991180 8 h1 FALSE
## 5 12.098570 12.233131 12 h2 FALSEin75 <- map(paramsgam, calc_HPDI, prob=0.75)## Joining, by = c("c.1_model", "c.2_model", "beta_model", "h1_model", "h2_model",
## "lp___model")## Warning: Column `c.1_model` has different attributes on LHS and RHS of join## Warning: Column `c.2_model` has different attributes on LHS and RHS of join## Warning: Column `beta_model` has different attributes on LHS and RHS of join## Warning: Column `h1_model` has different attributes on LHS and RHS of join## Warning: Column `h2_model` has different attributes on LHS and RHS of join## Warning: Column `lp___model` has different attributes on LHS and RHS of join## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")## Warning: Column `c.1_model` has different attributes on LHS and RHS of join## Warning: Column `c.2_model` has different attributes on LHS and RHS of join## Warning: Column `beta_model` has different attributes on LHS and RHS of join## Warning: Column `h1_model` has different attributes on LHS and RHS of join## Warning: Column `h2_model` has different attributes on LHS and RHS of join## Warning: Column `lp___model` has different attributes on LHS and RHS of join## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")## Warning: Column `c.1_model` has different attributes on LHS and RHS of join## Warning: Column `c.2_model` has different attributes on LHS and RHS of join## Warning: Column `beta_model` has different attributes on LHS and RHS of join## Warning: Column `h1_model` has different attributes on LHS and RHS of join## Warning: Column `h2_model` has different attributes on LHS and RHS of join## Warning: Column `lp___model` has different attributes on LHS and RHS of join## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")## Warning: Column `c.1_model` has different attributes on LHS and RHS of join## Warning: Column `c.2_model` has different attributes on LHS and RHS of join## Warning: Column `beta_model` has different attributes on LHS and RHS of join## Warning: Column `h1_model` has different attributes on LHS and RHS of join## Warning: Column `h2_model` has different attributes on LHS and RHS of join## Warning: Column `lp___model` has different attributes on LHS and RHS of join## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")## Warning: Column `c.1_model` has different attributes on LHS and RHS of join## Warning: Column `c.2_model` has different attributes on LHS and RHS of join## Warning: Column `beta_model` has different attributes on LHS and RHS of join## Warning: Column `h1_model` has different attributes on LHS and RHS of join## Warning: Column `h2_model` has different attributes on LHS and RHS of join## Warning: Column `lp___model` has different attributes on LHS and RHS of join## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")## Warning: Column `c.1_model` has different attributes on LHS and RHS of join## Warning: Column `c.2_model` has different attributes on LHS and RHS of join## Warning: Column `beta_model` has different attributes on LHS and RHS of join## Warning: Column `h1_model` has different attributes on LHS and RHS of join## Warning: Column `h2_model` has different attributes on LHS and RHS of join## Warning: Column `lp___model` has different attributes on LHS and RHS of join## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")## Warning: Column `c.1_model` has different attributes on LHS and RHS of join## Warning: Column `c.2_model` has different attributes on LHS and RHS of join## Warning: Column `beta_model` has different attributes on LHS and RHS of join## Warning: Column `h1_model` has different attributes on LHS and RHS of join## Warning: Column `h2_model` has different attributes on LHS and RHS of join## Warning: Column `lp___model` has different attributes on LHS and RHS of join## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")## Warning: Column `c.1_model` has different attributes on LHS and RHS of join## Warning: Column `c.2_model` has different attributes on LHS and RHS of join## Warning: Column `beta_model` has different attributes on LHS and RHS of join## Warning: Column `h1_model` has different attributes on LHS and RHS of join## Warning: Column `h2_model` has different attributes on LHS and RHS of join## Warning: Column `lp___model` has different attributes on LHS and RHS of join## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")## Warning: Column `c.1_model` has different attributes on LHS and RHS of join## Warning: Column `c.2_model` has different attributes on LHS and RHS of join## Warning: Column `beta_model` has different attributes on LHS and RHS of join## Warning: Column `h1_model` has different attributes on LHS and RHS of join## Warning: Column `h2_model` has different attributes on LHS and RHS of join## Warning: Column `lp___model` has different attributes on LHS and RHS of join## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")## Warning: Column `c.1_model` has different attributes on LHS and RHS of join## Warning: Column `c.2_model` has different attributes on LHS and RHS of join## Warning: Column `beta_model` has different attributes on LHS and RHS of join## Warning: Column `h1_model` has different attributes on LHS and RHS of join## Warning: Column `h2_model` has different attributes on LHS and RHS of join## Warning: Column `lp___model` has different attributes on LHS and RHS of join## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")## Warning: Column `c.1_model` has different attributes on LHS and RHS of join## Warning: Column `c.2_model` has different attributes on LHS and RHS of join## Warning: Column `beta_model` has different attributes on LHS and RHS of join## Warning: Column `h1_model` has different attributes on LHS and RHS of join## Warning: Column `h2_model` has different attributes on LHS and RHS of join## Warning: Column `lp___model` has different attributes on LHS and RHS of join## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")## Warning: Column `c.1_model` has different attributes on LHS and RHS of join## Warning: Column `c.2_model` has different attributes on LHS and RHS of join## Warning: Column `beta_model` has different attributes on LHS and RHS of join## Warning: Column `h1_model` has different attributes on LHS and RHS of join## Warning: Column `h2_model` has different attributes on LHS and RHS of join## Warning: Column `lp___model` has different attributes on LHS and RHS of join## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")## Warning: Column `c.1_model` has different attributes on LHS and RHS of join## Warning: Column `c.2_model` has different attributes on LHS and RHS of join## Warning: Column `beta_model` has different attributes on LHS and RHS of join## Warning: Column `h1_model` has different attributes on LHS and RHS of join## Warning: Column `h2_model` has different attributes on LHS and RHS of join## Warning: Column `lp___model` has different attributes on LHS and RHS of join## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")## Warning: Column `c.1_model` has different attributes on LHS and RHS of join## Warning: Column `c.2_model` has different attributes on LHS and RHS of join## Warning: Column `beta_model` has different attributes on LHS and RHS of join## Warning: Column `h1_model` has different attributes on LHS and RHS of join## Warning: Column `h2_model` has different attributes on LHS and RHS of join## Warning: Column `lp___model` has different attributes on LHS and RHS of join## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")## Warning: Column `c.1_model` has different attributes on LHS and RHS of join## Warning: Column `c.2_model` has different attributes on LHS and RHS of join## Warning: Column `beta_model` has different attributes on LHS and RHS of join## Warning: Column `h1_model` has different attributes on LHS and RHS of join## Warning: Column `h2_model` has different attributes on LHS and RHS of join## Warning: Column `lp___model` has different attributes on LHS and RHS of join## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")## Warning: Column `c.1_model` has different attributes on LHS and RHS of join## Warning: Column `c.2_model` has different attributes on LHS and RHS of join## Warning: Column `beta_model` has different attributes on LHS and RHS of join## Warning: Column `h1_model` has different attributes on LHS and RHS of join## Warning: Column `h2_model` has different attributes on LHS and RHS of join## Warning: Column `lp___model` has different attributes on LHS and RHS of join## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")## Warning: Column `c.1_model` has different attributes on LHS and RHS of join## Warning: Column `c.2_model` has different attributes on LHS and RHS of join## Warning: Column `beta_model` has different attributes on LHS and RHS of join## Warning: Column `h1_model` has different attributes on LHS and RHS of join## Warning: Column `h2_model` has different attributes on LHS and RHS of join## Warning: Column `lp___model` has different attributes on LHS and RHS of join## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")## Warning: Column `c.1_model` has different attributes on LHS and RHS of join## Warning: Column `c.2_model` has different attributes on LHS and RHS of join## Warning: Column `beta_model` has different attributes on LHS and RHS of join## Warning: Column `h1_model` has different attributes on LHS and RHS of join## Warning: Column `h2_model` has different attributes on LHS and RHS of join## Warning: Column `lp___model` has different attributes on LHS and RHS of join## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")print(in75)## $m1
## low high true params inint
## 1 3.2501236 3.9272662 4.0 c.1 FALSE
## 2 4.9079239 5.6276046 6.0 c.2 FALSE
## 3 0.4066015 0.4748422 0.5 beta FALSE
## 4 7.8773879 8.2708234 8.0 h1 TRUE
## 5 11.7804592 12.1964284 12.0 h2 TRUE
##
## $m2
## low high true params inint
## 1 3.1747532 3.8290328 4.0 c.1 FALSE
## 2 4.9277306 5.6227433 6.0 c.2 FALSE
## 3 0.4031203 0.4692065 0.5 beta FALSE
## 4 7.8930745 8.2767062 8.0 h1 TRUE
## 5 11.8244379 12.2359328 12.0 h2 TRUE
##
## $m3
## low high true params inint
## 1 3.1904716 3.8115469 4.0 c.1 FALSE
## 2 4.9024314 5.5857277 6.0 c.2 FALSE
## 3 0.4086681 0.4723668 0.5 beta FALSE
## 4 7.8722504 8.2621884 8.0 h1 TRUE
## 5 11.8230640 12.2186223 12.0 h2 TRUE
##
## $m4
## low high true params inint
## 1 3.4753182 4.1142259 4.0 c.1 TRUE
## 2 5.2798238 5.9367518 6.0 c.2 FALSE
## 3 0.4363857 0.4991233 0.5 beta FALSE
## 4 8.0089769 8.3662952 8.0 h1 FALSE
## 5 11.7793695 12.1459087 12.0 h2 TRUE
##
## $m5
## low high true params inint
## 1 3.4433367 4.0945892 4.0 c.1 TRUE
## 2 5.1325837 5.8360014 6.0 c.2 FALSE
## 3 0.4246286 0.4904984 0.5 beta FALSE
## 4 8.0189219 8.3836955 8.0 h1 FALSE
## 5 11.7668734 12.1422131 12.0 h2 TRUE
##
## $m6
## low high true params inint
## 1 3.370690 4.0211574 4.0 c.1 TRUE
## 2 5.186302 5.8966790 6.0 c.2 FALSE
## 3 0.418465 0.4850303 0.5 beta FALSE
## 4 7.999701 8.3563466 8.0 h1 TRUE
## 5 11.791000 12.1658031 12.0 h2 TRUE
##
## $m7
## low high true params inint
## 1 3.9334486 4.5523476 4.0 c.1 TRUE
## 2 5.7414360 6.4048148 6.0 c.2 TRUE
## 3 0.4744432 0.5376042 0.5 beta TRUE
## 4 8.2128315 8.5227315 8.0 h1 FALSE
## 5 11.6305653 11.9513433 12.0 h2 FALSE
##
## $m8
## low high true params inint
## 1 3.9378608 4.557596 4.0 c.1 TRUE
## 2 5.7729105 6.447308 6.0 c.2 TRUE
## 3 0.4702384 0.533441 0.5 beta TRUE
## 4 8.2315496 8.543833 8.0 h1 FALSE
## 5 11.6784857 11.993660 12.0 h2 FALSE
##
## $m9
## low high true params inint
## 1 3.9731254 4.5753073 4.0 c.1 TRUE
## 2 5.6989545 6.3508008 6.0 c.2 TRUE
## 3 0.4805663 0.5419001 0.5 beta TRUE
## 4 8.1834663 8.5009195 8.0 h1 FALSE
## 5 11.6774571 11.9984272 12.0 h2 FALSE
##
## $m10
## low high true params inint
## 1 11.159384 12.494009 16 c.1 FALSE
## 2 18.014987 19.826380 24 c.2 FALSE
## 3 1.448803 1.602843 2 beta FALSE
## 4 7.672810 7.841917 8 h1 FALSE
## 5 12.192592 12.356752 12 h2 FALSE
##
## $m11
## low high true params inint
## 1 11.185234 12.402412 16 c.1 FALSE
## 2 17.852498 19.520479 24 c.2 FALSE
## 3 1.445014 1.587621 2 beta FALSE
## 4 7.656411 7.817452 8 h1 FALSE
## 5 12.170768 12.333943 12 h2 FALSE
##
## $m12
## low high true params inint
## 1 10.848001 12.049589 16 c.1 FALSE
## 2 17.692949 19.370861 24 c.2 FALSE
## 3 1.421928 1.564056 2 beta FALSE
## 4 7.664332 7.829884 8 h1 FALSE
## 5 12.197229 12.360644 12 h2 FALSE
##
## $m13
## low high true params inint
## 1 12.524896 13.899434 16 c.1 FALSE
## 2 19.839749 21.731258 24 c.2 FALSE
## 3 1.604186 1.765372 2 beta FALSE
## 4 7.760157 7.913378 8 h1 FALSE
## 5 12.146768 12.297673 12 h2 FALSE
##
## $m14
## low high true params inint
## 1 12.320277 13.634744 16 c.1 FALSE
## 2 19.447600 21.265906 24 c.2 FALSE
## 3 1.592442 1.745851 2 beta FALSE
## 4 7.735040 7.883800 8 h1 FALSE
## 5 12.148636 12.299303 12 h2 FALSE
##
## $m15
## low high true params inint
## 1 12.213129 13.531782 16 c.1 FALSE
## 2 19.267255 21.137336 24 c.2 FALSE
## 3 1.575907 1.733495 2 beta FALSE
## 4 7.729892 7.877785 8 h1 FALSE
## 5 12.143976 12.292382 12 h2 FALSE
##
## $m16
## low high true params inint
## 1 15.299240 16.740585 16 c.1 TRUE
## 2 23.408511 25.483513 24 c.2 TRUE
## 3 1.930559 2.106695 2 beta TRUE
## 4 7.854474 7.987955 8 h1 FALSE
## 5 12.089570 12.224809 12 h2 FALSE
##
## $m17
## low high true params inint
## 1 15.322033 16.795263 16 c.1 TRUE
## 2 23.497265 25.605068 24 c.2 TRUE
## 3 1.939864 2.115915 2 beta TRUE
## 4 7.845753 7.980946 8 h1 FALSE
## 5 12.070039 12.206998 12 h2 FALSE
##
## $m18
## low high true params inint
## 1 14.975978 16.394557 16 c.1 TRUE
## 2 23.018115 24.980390 24 c.2 TRUE
## 3 1.907683 2.074526 2 beta TRUE
## 4 7.856500 7.991180 8 h1 FALSE
## 5 12.098570 12.233131 12 h2 FALSEin90 <- map(paramsgam, calc_HPDI, prob=0.90)## Joining, by = c("c.1_model", "c.2_model", "beta_model", "h1_model", "h2_model",
## "lp___model")## Warning: Column `c.1_model` has different attributes on LHS and RHS of join## Warning: Column `c.2_model` has different attributes on LHS and RHS of join## Warning: Column `beta_model` has different attributes on LHS and RHS of join## Warning: Column `h1_model` has different attributes on LHS and RHS of join## Warning: Column `h2_model` has different attributes on LHS and RHS of join## Warning: Column `lp___model` has different attributes on LHS and RHS of join## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")## Warning: Column `c.1_model` has different attributes on LHS and RHS of join## Warning: Column `c.2_model` has different attributes on LHS and RHS of join## Warning: Column `beta_model` has different attributes on LHS and RHS of join## Warning: Column `h1_model` has different attributes on LHS and RHS of join## Warning: Column `h2_model` has different attributes on LHS and RHS of join## Warning: Column `lp___model` has different attributes on LHS and RHS of join## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")## Warning: Column `c.1_model` has different attributes on LHS and RHS of join## Warning: Column `c.2_model` has different attributes on LHS and RHS of join## Warning: Column `beta_model` has different attributes on LHS and RHS of join## Warning: Column `h1_model` has different attributes on LHS and RHS of join## Warning: Column `h2_model` has different attributes on LHS and RHS of join## Warning: Column `lp___model` has different attributes on LHS and RHS of join## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")## Warning: Column `c.1_model` has different attributes on LHS and RHS of join## Warning: Column `c.2_model` has different attributes on LHS and RHS of join## Warning: Column `beta_model` has different attributes on LHS and RHS of join## Warning: Column `h1_model` has different attributes on LHS and RHS of join## Warning: Column `h2_model` has different attributes on LHS and RHS of join## Warning: Column `lp___model` has different attributes on LHS and RHS of join## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")## Warning: Column `c.1_model` has different attributes on LHS and RHS of join## Warning: Column `c.2_model` has different attributes on LHS and RHS of join## Warning: Column `beta_model` has different attributes on LHS and RHS of join## Warning: Column `h1_model` has different attributes on LHS and RHS of join## Warning: Column `h2_model` has different attributes on LHS and RHS of join## Warning: Column `lp___model` has different attributes on LHS and RHS of join## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")## Warning: Column `c.1_model` has different attributes on LHS and RHS of join## Warning: Column `c.2_model` has different attributes on LHS and RHS of join## Warning: Column `beta_model` has different attributes on LHS and RHS of join## Warning: Column `h1_model` has different attributes on LHS and RHS of join## Warning: Column `h2_model` has different attributes on LHS and RHS of join## Warning: Column `lp___model` has different attributes on LHS and RHS of join## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")## Warning: Column `c.1_model` has different attributes on LHS and RHS of join## Warning: Column `c.2_model` has different attributes on LHS and RHS of join## Warning: Column `beta_model` has different attributes on LHS and RHS of join## Warning: Column `h1_model` has different attributes on LHS and RHS of join## Warning: Column `h2_model` has different attributes on LHS and RHS of join## Warning: Column `lp___model` has different attributes on LHS and RHS of join## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")## Warning: Column `c.1_model` has different attributes on LHS and RHS of join## Warning: Column `c.2_model` has different attributes on LHS and RHS of join## Warning: Column `beta_model` has different attributes on LHS and RHS of join## Warning: Column `h1_model` has different attributes on LHS and RHS of join## Warning: Column `h2_model` has different attributes on LHS and RHS of join## Warning: Column `lp___model` has different attributes on LHS and RHS of join## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")## Warning: Column `c.1_model` has different attributes on LHS and RHS of join## Warning: Column `c.2_model` has different attributes on LHS and RHS of join## Warning: Column `beta_model` has different attributes on LHS and RHS of join## Warning: Column `h1_model` has different attributes on LHS and RHS of join## Warning: Column `h2_model` has different attributes on LHS and RHS of join## Warning: Column `lp___model` has different attributes on LHS and RHS of join## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")## Warning: Column `c.1_model` has different attributes on LHS and RHS of join## Warning: Column `c.2_model` has different attributes on LHS and RHS of join## Warning: Column `beta_model` has different attributes on LHS and RHS of join## Warning: Column `h1_model` has different attributes on LHS and RHS of join## Warning: Column `h2_model` has different attributes on LHS and RHS of join## Warning: Column `lp___model` has different attributes on LHS and RHS of join## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")## Warning: Column `c.1_model` has different attributes on LHS and RHS of join## Warning: Column `c.2_model` has different attributes on LHS and RHS of join## Warning: Column `beta_model` has different attributes on LHS and RHS of join## Warning: Column `h1_model` has different attributes on LHS and RHS of join## Warning: Column `h2_model` has different attributes on LHS and RHS of join## Warning: Column `lp___model` has different attributes on LHS and RHS of join## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")## Warning: Column `c.1_model` has different attributes on LHS and RHS of join## Warning: Column `c.2_model` has different attributes on LHS and RHS of join## Warning: Column `beta_model` has different attributes on LHS and RHS of join## Warning: Column `h1_model` has different attributes on LHS and RHS of join## Warning: Column `h2_model` has different attributes on LHS and RHS of join## Warning: Column `lp___model` has different attributes on LHS and RHS of join## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")## Warning: Column `c.1_model` has different attributes on LHS and RHS of join## Warning: Column `c.2_model` has different attributes on LHS and RHS of join## Warning: Column `beta_model` has different attributes on LHS and RHS of join## Warning: Column `h1_model` has different attributes on LHS and RHS of join## Warning: Column `h2_model` has different attributes on LHS and RHS of join## Warning: Column `lp___model` has different attributes on LHS and RHS of join## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")## Warning: Column `c.1_model` has different attributes on LHS and RHS of join## Warning: Column `c.2_model` has different attributes on LHS and RHS of join## Warning: Column `beta_model` has different attributes on LHS and RHS of join## Warning: Column `h1_model` has different attributes on LHS and RHS of join## Warning: Column `h2_model` has different attributes on LHS and RHS of join## Warning: Column `lp___model` has different attributes on LHS and RHS of join## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")## Warning: Column `c.1_model` has different attributes on LHS and RHS of join## Warning: Column `c.2_model` has different attributes on LHS and RHS of join## Warning: Column `beta_model` has different attributes on LHS and RHS of join## Warning: Column `h1_model` has different attributes on LHS and RHS of join## Warning: Column `h2_model` has different attributes on LHS and RHS of join## Warning: Column `lp___model` has different attributes on LHS and RHS of join## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")## Warning: Column `c.1_model` has different attributes on LHS and RHS of join## Warning: Column `c.2_model` has different attributes on LHS and RHS of join## Warning: Column `beta_model` has different attributes on LHS and RHS of join## Warning: Column `h1_model` has different attributes on LHS and RHS of join## Warning: Column `h2_model` has different attributes on LHS and RHS of join## Warning: Column `lp___model` has different attributes on LHS and RHS of join## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")## Warning: Column `c.1_model` has different attributes on LHS and RHS of join## Warning: Column `c.2_model` has different attributes on LHS and RHS of join## Warning: Column `beta_model` has different attributes on LHS and RHS of join## Warning: Column `h1_model` has different attributes on LHS and RHS of join## Warning: Column `h2_model` has different attributes on LHS and RHS of join## Warning: Column `lp___model` has different attributes on LHS and RHS of join## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")## Warning: Column `c.1_model` has different attributes on LHS and RHS of join## Warning: Column `c.2_model` has different attributes on LHS and RHS of join## Warning: Column `beta_model` has different attributes on LHS and RHS of join## Warning: Column `h1_model` has different attributes on LHS and RHS of join## Warning: Column `h2_model` has different attributes on LHS and RHS of join## Warning: Column `lp___model` has different attributes on LHS and RHS of join## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")print(in90)## $m1
## low high true params inint
## 1 3.2501236 3.9272662 4.0 c.1 FALSE
## 2 4.9079239 5.6276046 6.0 c.2 FALSE
## 3 0.4066015 0.4748422 0.5 beta FALSE
## 4 7.8773879 8.2708234 8.0 h1 TRUE
## 5 11.7804592 12.1964284 12.0 h2 TRUE
##
## $m2
## low high true params inint
## 1 3.1747532 3.8290328 4.0 c.1 FALSE
## 2 4.9277306 5.6227433 6.0 c.2 FALSE
## 3 0.4031203 0.4692065 0.5 beta FALSE
## 4 7.8930745 8.2767062 8.0 h1 TRUE
## 5 11.8244379 12.2359328 12.0 h2 TRUE
##
## $m3
## low high true params inint
## 1 3.1904716 3.8115469 4.0 c.1 FALSE
## 2 4.9024314 5.5857277 6.0 c.2 FALSE
## 3 0.4086681 0.4723668 0.5 beta FALSE
## 4 7.8722504 8.2621884 8.0 h1 TRUE
## 5 11.8230640 12.2186223 12.0 h2 TRUE
##
## $m4
## low high true params inint
## 1 3.4753182 4.1142259 4.0 c.1 TRUE
## 2 5.2798238 5.9367518 6.0 c.2 FALSE
## 3 0.4363857 0.4991233 0.5 beta FALSE
## 4 8.0089769 8.3662952 8.0 h1 FALSE
## 5 11.7793695 12.1459087 12.0 h2 TRUE
##
## $m5
## low high true params inint
## 1 3.4433367 4.0945892 4.0 c.1 TRUE
## 2 5.1325837 5.8360014 6.0 c.2 FALSE
## 3 0.4246286 0.4904984 0.5 beta FALSE
## 4 8.0189219 8.3836955 8.0 h1 FALSE
## 5 11.7668734 12.1422131 12.0 h2 TRUE
##
## $m6
## low high true params inint
## 1 3.370690 4.0211574 4.0 c.1 TRUE
## 2 5.186302 5.8966790 6.0 c.2 FALSE
## 3 0.418465 0.4850303 0.5 beta FALSE
## 4 7.999701 8.3563466 8.0 h1 TRUE
## 5 11.791000 12.1658031 12.0 h2 TRUE
##
## $m7
## low high true params inint
## 1 3.9334486 4.5523476 4.0 c.1 TRUE
## 2 5.7414360 6.4048148 6.0 c.2 TRUE
## 3 0.4744432 0.5376042 0.5 beta TRUE
## 4 8.2128315 8.5227315 8.0 h1 FALSE
## 5 11.6305653 11.9513433 12.0 h2 FALSE
##
## $m8
## low high true params inint
## 1 3.9378608 4.557596 4.0 c.1 TRUE
## 2 5.7729105 6.447308 6.0 c.2 TRUE
## 3 0.4702384 0.533441 0.5 beta TRUE
## 4 8.2315496 8.543833 8.0 h1 FALSE
## 5 11.6784857 11.993660 12.0 h2 FALSE
##
## $m9
## low high true params inint
## 1 3.9731254 4.5753073 4.0 c.1 TRUE
## 2 5.6989545 6.3508008 6.0 c.2 TRUE
## 3 0.4805663 0.5419001 0.5 beta TRUE
## 4 8.1834663 8.5009195 8.0 h1 FALSE
## 5 11.6774571 11.9984272 12.0 h2 FALSE
##
## $m10
## low high true params inint
## 1 11.159384 12.494009 16 c.1 FALSE
## 2 18.014987 19.826380 24 c.2 FALSE
## 3 1.448803 1.602843 2 beta FALSE
## 4 7.672810 7.841917 8 h1 FALSE
## 5 12.192592 12.356752 12 h2 FALSE
##
## $m11
## low high true params inint
## 1 11.185234 12.402412 16 c.1 FALSE
## 2 17.852498 19.520479 24 c.2 FALSE
## 3 1.445014 1.587621 2 beta FALSE
## 4 7.656411 7.817452 8 h1 FALSE
## 5 12.170768 12.333943 12 h2 FALSE
##
## $m12
## low high true params inint
## 1 10.848001 12.049589 16 c.1 FALSE
## 2 17.692949 19.370861 24 c.2 FALSE
## 3 1.421928 1.564056 2 beta FALSE
## 4 7.664332 7.829884 8 h1 FALSE
## 5 12.197229 12.360644 12 h2 FALSE
##
## $m13
## low high true params inint
## 1 12.524896 13.899434 16 c.1 FALSE
## 2 19.839749 21.731258 24 c.2 FALSE
## 3 1.604186 1.765372 2 beta FALSE
## 4 7.760157 7.913378 8 h1 FALSE
## 5 12.146768 12.297673 12 h2 FALSE
##
## $m14
## low high true params inint
## 1 12.320277 13.634744 16 c.1 FALSE
## 2 19.447600 21.265906 24 c.2 FALSE
## 3 1.592442 1.745851 2 beta FALSE
## 4 7.735040 7.883800 8 h1 FALSE
## 5 12.148636 12.299303 12 h2 FALSE
##
## $m15
## low high true params inint
## 1 12.213129 13.531782 16 c.1 FALSE
## 2 19.267255 21.137336 24 c.2 FALSE
## 3 1.575907 1.733495 2 beta FALSE
## 4 7.729892 7.877785 8 h1 FALSE
## 5 12.143976 12.292382 12 h2 FALSE
##
## $m16
## low high true params inint
## 1 15.299240 16.740585 16 c.1 TRUE
## 2 23.408511 25.483513 24 c.2 TRUE
## 3 1.930559 2.106695 2 beta TRUE
## 4 7.854474 7.987955 8 h1 FALSE
## 5 12.089570 12.224809 12 h2 FALSE
##
## $m17
## low high true params inint
## 1 15.322033 16.795263 16 c.1 TRUE
## 2 23.497265 25.605068 24 c.2 TRUE
## 3 1.939864 2.115915 2 beta TRUE
## 4 7.845753 7.980946 8 h1 FALSE
## 5 12.070039 12.206998 12 h2 FALSE
##
## $m18
## low high true params inint
## 1 14.975978 16.394557 16 c.1 TRUE
## 2 23.018115 24.980390 24 c.2 TRUE
## 3 1.907683 2.074526 2 beta TRUE
## 4 7.856500 7.991180 8 h1 FALSE
## 5 12.098570 12.233131 12 h2 FALSEGamma prior recaptures parameters relatively well here and there are no obvious model fitting problems. However, it pretty much always overestimates beta and cutpoints (while accurately identifying the inflection points). It’s not especially sensitive to the prior on \(\beta\) either. An \(\mathrm{exponential}\) prior with rates of 1,2,3, and 6 on \(\beta\) performed similarly.